{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# relax_dynamic calculation style\n",
    "\n",
    "**Lucas M. Hale**, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), *Materials Science and Engineering Division, NIST*.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The relax_dynamic calculation style dynamically relaxes an atomic configuration for a specified number of timesteps.  Upon completion, the mean, $\\langle X \\rangle$, and standard deviation, $\\sigma_X$, of all thermo properties, $X$, are computed for a specified range of times.  This method is meant to measure equilibrium properties of bulk materials, both at zero K and at various temperatures.\n",
    "\n",
    "### Version notes\n",
    "\n",
    "- 2018-07-09: Notebook added.\n",
    "- 2019-07-30: Description updated and small changes due to iprPy version.\n",
    "- 2020-05-22: Version 0.10 update - potentials now loaded from database.\n",
    "- 2020-09-22: Setup and parameter definition streamlined.\n",
    "\n",
    "### Additional dependencies\n",
    "\n",
    "### Disclaimers\n",
    "\n",
    "- [NIST disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)\n",
    "- The calculation reports the standard deviation, $\\sigma_X$ of the measured properties not the standard error of the mean, $\\sigma_{\\langle X \\rangle}$.  The two are related to each other according to $\\sigma_{\\langle X \\rangle} = \\sigma_X \\sqrt{\\frac{C}{N}}$, where $N$ is the number of samples taken of $X$, and $C$ is a statistical inefficiency due to the autocorrelation of the measurements with time.  Obtaining a proper estimate of $\\sigma_{\\langle X \\rangle}$ requires either estimating $C$ from the raw thermo data (not done here), or only taking measurements sporadically to ensure the samples are independent.\n",
    "- Good (low error) results requires running large simulations for a long time.  The reasons for this are:\n",
    "  - Systems have to be large enough to avoid issues with fluctuations across the periodic boundaries.\n",
    "  - Runs must first let the systems equilibrate before meaningful measurements can be taken.\n",
    "  - The standard deviation, $\\sigma$, of thermo properties is proportional to the number of atoms, $N_a$ as $\\sigma \\propto \\frac{1}{\\sqrt{N_a}}$.\n",
    "  - The standard error, $\\sigma_x$ of thermo properties is proportional to the number of samples taken, $N$ as $\\sigma_x \\propto \\frac{1}{\\sqrt{N}}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method and Theory\n",
    "\n",
    "An initial system (and corresponding unit cell system) is supplied with box dimensions, $a_i^0$, close to the equilibrium values. A LAMMPS simulation then integrates the atomic positions and velocities for a specified number of timesteps.\n",
    "\n",
    "The calculation script allows for the use of different integration methods:\n",
    "\n",
    "- nve integrates atomic positions without changing box dimensions or the system's total energy.\n",
    "\n",
    "- npt integrates atomic positions and applies Nose-Hoover style thermostat and barostat (equilibriate to specified T and P).\n",
    "\n",
    "- nvt integrates atomic positions and applies Nose-Hoover style thermostat (equilibriate to specified T).\n",
    "\n",
    "- nph integrates atomic positions and applies Nose-Hoover style barostat (equilibriate to specified P).\n",
    "\n",
    "- nve+l integrates atomic positions and applies Langevin style thermostat (equilibriate to specified T).\n",
    "\n",
    "- nph+l integrates atomic positions and applies Nose-Hoover style barostat and Langevin style thermostat (equilibriate to specified T and P).\n",
    "\n",
    "__Notes__ on the different control schemes:\n",
    "\n",
    "- The Nose-Hoover barostat works by rescaling the box dimensions according to the measured system pressures.\n",
    "\n",
    "- The Nose-Hoover thermostat works by rescaling the atomic velocities according to the measured system temperature (kinetic energy). Cannot be used with a temperature of 0 K.\n",
    "\n",
    "- The Langevin thermostat works by modifying the forces on all atoms with both a dampener and a random temperature dependent fluctuation. Used at 0 K, only the force dampener is applied.\n",
    "\n",
    "__Notes__ on run parameter values. The proper time to reach equilibrium (equilsteps), and sample frequency to ensure uncorrelated measurements (thermosteps) is simulation dependent. They can be influenced by the potential, timestep size, crystal structure, integration method, presence of defects, etc. The default values of equilsteps = 20,000 and thermosteps = 100 are based on general rule-of-thumb estimates for bulk crystals and EAM potentials, and may or may not be adequate.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1. Library imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import libraries needed by the calculation. The external libraries used are:\n",
    "\n",
    "- [numpy](http://www.numpy.org/)\n",
    "\n",
    "- [atomman](https://github.com/usnistgov/atomman)\n",
    "\n",
    "- [iprPy](https://github.com/usnistgov/iprPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook last executed on 2020-09-22 using iprPy version 0.10.2\n"
     ]
    }
   ],
   "source": [
    "# Standard library imports\n",
    "from pathlib import Path\n",
    "import os\n",
    "import datetime\n",
    "import random\n",
    "from copy import deepcopy\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np\n",
    "\n",
    "# https://github.com/usnistgov/atomman \n",
    "import atomman as am\n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# https://github.com/usnistgov/iprPy\n",
    "import iprPy\n",
    "\n",
    "print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. Default calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify calculation style\n",
    "calc_style = 'relax_dynamic'\n",
    "\n",
    "# If workingdir is already set, then do nothing (already in correct folder)\n",
    "try:\n",
    "    workingdir = workingdir\n",
    "\n",
    "# Change to workingdir if not already there\n",
    "except:\n",
    "    workingdir = Path('calculationfiles', calc_style)\n",
    "    if not workingdir.is_dir():\n",
    "        workingdir.mkdir(parents=True)\n",
    "    os.chdir(workingdir)\n",
    "    \n",
    "# Initialize connection to library\n",
    "library = iprPy.Library(load=['lammps_potentials'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign values for the calculation's run parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1. Specify system-specific paths\n",
    "\n",
    "- __lammps_command__ is the LAMMPS command to use (required).\n",
    "\n",
    "- __mpi_command__ MPI command for running LAMMPS in parallel. A value of None will run simulations serially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lammps_command = 'lmp_mpi'\n",
    "#mpi_command = None\n",
    "mpi_command = 'C:/Program Files/MPICH2/bin/mpiexec.exe -localonly 6'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2. Load interatomic potential\n",
    "\n",
    "- __potential_name__ gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.  \n",
    "\n",
    "- __potential__ is an atomman.lammps.Potential object (required)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'\n",
    "\n",
    "# Retrieve potential and parameter file(s)\n",
    "potential = library.get_lammps_potential(id=potential_name, getfiles=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3. Load initial unit cell system\n",
    "\n",
    "- __ucell__ is an atomman.System representing a fundamental unit cell of the system (required).  Here, this is generated using the load parameters and symbols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.500,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  3.500,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  3.500]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.750 |   1.750\n",
      "      2 |       1 |   1.750 |   0.000 |   1.750\n",
      "      3 |       1 |   1.750 |   1.750 |   0.000\n"
     ]
    }
   ],
   "source": [
    "# Create ucell by loading prototype record\n",
    "ucell = am.load('prototype', 'A1--Cu--fcc', symbols='Ni', a=3.5)\n",
    "\n",
    "print(ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.4. Modify system\n",
    "\n",
    "- __sizemults__ list of three integers specifying how many times the ucell vectors of $a$, $b$ and $c$ are replicated in creating system.\n",
    "\n",
    "- __system__ is an atomman.System to perform the scan on (required). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# of atoms in system = 4000\n"
     ]
    }
   ],
   "source": [
    "sizemults = [10, 10, 10]\n",
    "\n",
    "# Generate system by supersizing ucell\n",
    "system = ucell.supersize(*sizemults)\n",
    "print('# of atoms in system =', system.natoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.5. Specify calculation-specific run parameters\n",
    "\n",
    "- __pressure_xx__ gives the xx component of the pressure to equilibriate the system to (npt, nph, and nph+l styles).\n",
    "\n",
    "- __pressure_yy__ gives the yy component of the pressure to equilibriate the system to (npt, nph, and nph+l styles).\n",
    "\n",
    "- __pressure_zz__ gives the zz component of the pressure to equilibriate the system to (npt, nph, and nph+l styles).\n",
    "\n",
    "- __pressure_xy__ gives the xy component of the pressure to equilibriate the system to (npt, nph, and nph+l styles).\n",
    "\n",
    "- __pressure_xz__ gives the xz component of the pressure to equilibriate the system to (npt, nph, and nph+l styles).\n",
    "\n",
    "- __pressure_yz__ gives the yz component of the pressure to equilibriate the system to (npt, nph, and nph+l styles).\n",
    "\n",
    "- __temperature__ gives the temperature to equilibriate the system to (nvt, npt, nve+l, and nph+l styles).\n",
    "\n",
    "- __integrator__ specifies the integrator style to use. Default value is 'nph+l' for temperature = 0, and 'npt' otherwise.\n",
    "\n",
    "- __runsteps__ is the total number of integration timesteps to perform. Default value is 220000.\n",
    "\n",
    "- __thermosteps__ specifies to output thermo values every this many timesteps. Default value is 100.\n",
    "    \n",
    "- __dumpsteps__ specifies to output dump files every this many timesteps. Default value is runsteps (only first and last steps are outputted as dump files).\n",
    "    \n",
    "- __equilsteps__ is the number of timesteps to equilibriate the system for. Only thermo values associated with timesteps greater than equilsteps will be included in the mean and standard deviation calculations. Default value is 20000. \n",
    "\n",
    "- __randomseed__ specifies a random number seed used to generate the initial atomic velocities and the Langevin thermostat fluctuations. Default value generates a new random integer every time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "pressure_xx = uc.set_in_units(0.0, 'GPa')\n",
    "pressure_yy = uc.set_in_units(0.0, 'GPa')\n",
    "pressure_zz = uc.set_in_units(0.0, 'GPa')\n",
    "pressure_xy = uc.set_in_units(0.0, 'GPa')\n",
    "pressure_xz = uc.set_in_units(0.0, 'GPa')\n",
    "pressure_yz = uc.set_in_units(0.0, 'GPa')\n",
    "temperature = 300.0\n",
    "integrator = 'npt'\n",
    "runsteps = 220000\n",
    "thermosteps = 100\n",
    "dumpsteps = runsteps\n",
    "equilsteps = 20000\n",
    "randomseed = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define calculation function(s) and generate template LAMMPS script(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1. full_relax.template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('full_relax.template', 'w') as f:\n",
    "    f.write(\"\"\"#LAMMPS input script that performs a simple dynamic integration\n",
    "\n",
    "box tilt large\n",
    "\n",
    "<atomman_system_pair_info>\n",
    "\n",
    "change_box all triclinic\n",
    "\n",
    "compute pe all pe/atom\n",
    "compute ke all ke/atom\n",
    "compute stress all stress/atom <stressterm>\n",
    "\n",
    "thermo <thermosteps>\n",
    "thermo_style custom step temp pe ke etotal lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy\n",
    "thermo_modify format float %.13e\n",
    "timestep 0.001\n",
    "\n",
    "<integrator_info>\n",
    "\n",
    "dump dumpit all custom <dumpsteps> *.dump <dump_keys>\n",
    "dump_modify dumpit format <dump_modify_format>\n",
    "restart <runsteps> *.restart\n",
    "\n",
    "run <runsteps> upto\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2. integrator_info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrator_info(integrator=None, p_xx=0.0, p_yy=0.0, p_zz=0.0, p_xy=0.0,\n",
    "                    p_xz=0.0, p_yz=0.0, temperature=0.0, randomseed=None,\n",
    "                    units='metal'):\n",
    "    \"\"\"\n",
    "    Generates LAMMPS commands for velocity creation and fix integrators. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    integrator : str or None, optional\n",
    "        The integration method to use. Options are 'npt', 'nvt', 'nph',\n",
    "        'nve', 'nve+l', 'nph+l'. The +l options use Langevin thermostat.\n",
    "        (Default is None, which will use 'nph+l' for temperature == 0, and\n",
    "        'npt' otherwise.)\n",
    "    p_xx : float, optional\n",
    "        The value to relax the x tensile pressure component to (default is\n",
    "        0.0).\n",
    "    p_yy : float, optional\n",
    "        The value to relax the y tensile pressure component to (default is\n",
    "        0.0).\n",
    "    p_zz : float, optional\n",
    "        The value to relax the z tensile pressure component to (default is\n",
    "        0.0).\n",
    "    p_xy : float, optional\n",
    "        The value to relax the xy shear pressure component to (default is\n",
    "        0.0).\n",
    "    p_xz : float, optional\n",
    "        The value to relax the xz shear pressure component to (default is\n",
    "        0.0).\n",
    "    p_yz : float, optional\n",
    "        The value to relax the yz shear pressure component to (default is\n",
    "        0.0).\n",
    "    temperature : float, optional\n",
    "        The temperature to relax at (default is 0.0).\n",
    "    randomseed : int or None, optional\n",
    "        Random number seed used by LAMMPS in creating velocities and with\n",
    "        the Langevin thermostat.  (Default is None which will select a\n",
    "        random int between 1 and 900000000.)\n",
    "    units : str, optional\n",
    "        The LAMMPS units style to use (default is 'metal').\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    str\n",
    "        The generated LAMMPS input lines for velocity create and fix\n",
    "        integration commands.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Get lammps units\n",
    "    lammps_units = lmp.style.unit(units)\n",
    "    Px = uc.get_in_units(p_xx, lammps_units['pressure'])\n",
    "    Py = uc.get_in_units(p_yy, lammps_units['pressure'])\n",
    "    Pz = uc.get_in_units(p_zz, lammps_units['pressure'])\n",
    "    Pxy = uc.get_in_units(p_xy, lammps_units['pressure'])\n",
    "    Pxz = uc.get_in_units(p_xz, lammps_units['pressure'])\n",
    "    Pyz = uc.get_in_units(p_yz, lammps_units['pressure'])\n",
    "    T = temperature\n",
    "    \n",
    "    # Check temperature and set default integrator\n",
    "    if temperature == 0.0:\n",
    "        if integrator is None: integrator = 'nph+l'\n",
    "        assert integrator not in ['npt', 'nvt'], 'npt and nvt cannot run at 0 K'\n",
    "    elif temperature > 0:\n",
    "        if integrator is None: integrator = 'npt'\n",
    "    else:\n",
    "        raise ValueError('Temperature must be positive')\n",
    "    \n",
    "    # Set default randomseed\n",
    "    if randomseed is None: randomseed = random.randint(1, 900000000)\n",
    "    \n",
    "    if integrator == 'npt':\n",
    "        start_temp = T*2.+1\n",
    "        Tdamp = 100 * lmp.style.timestep(units)\n",
    "        Pdamp = 1000 * lmp.style.timestep(units)\n",
    "        int_info = '\\n'.join([\n",
    "                'velocity all create %f %i' % (start_temp, randomseed),\n",
    "                'fix npt all npt temp %f %f %f &' % (T, T, Tdamp),\n",
    "                '                x %f %f %f &' % (Px, Px, Pdamp),\n",
    "                '                y %f %f %f &' % (Py, Py, Pdamp),\n",
    "                '                z %f %f %f &' % (Pz, Pz, Pdamp),\n",
    "                '                xy %f %f %f &' % (Pxy, Pxy, Pdamp),\n",
    "                '                xz %f %f %f &' % (Pxz, Pxz, Pdamp),\n",
    "                '                yz %f %f %f' % (Pyz, Pyz, Pdamp),\n",
    "                ])\n",
    "    \n",
    "    elif integrator == 'nvt':\n",
    "        start_temp = T*2.+1\n",
    "        Tdamp = 100 * lmp.style.timestep(units)\n",
    "        int_info = '\\n'.join([\n",
    "                'velocity all create %f %i' % (start_temp, randomseed),\n",
    "                'fix nvt all nvt temp %f %f %f' % (T, T, Tdamp),\n",
    "                ])\n",
    "    \n",
    "    elif integrator == 'nph':\n",
    "        Pdamp = 1000 * lmp.style.timestep(units)\n",
    "        int_info = '\\n'.join([\n",
    "                'fix nph all nph x %f %f %f &' % (Px, Px, Pdamp),\n",
    "                '                y %f %f %f &' % (Py, Py, Pdamp),\n",
    "                '                z %f %f %f &' % (Pz, Pz, Pdamp),\n",
    "                '                xy %f %f %f &' % (Pxy, Pxy, Pdamp),\n",
    "                '                xz %f %f %f &' % (Pxz, Pxz, Pdamp),\n",
    "                '                yz %f %f %f' % (Pyz, Pyz, Pdamp),\n",
    "                ])\n",
    "    \n",
    "    elif integrator == 'nve':\n",
    "        int_info = 'fix nve all nve'\n",
    "        \n",
    "    elif integrator == 'nve+l':\n",
    "        start_temp = T*2.+1\n",
    "        Tdamp = 100 * lmp.style.timestep(units)\n",
    "        int_info = '\\n'.join([\n",
    "                'velocity all create %f %i' % (start_temp, randomseed),\n",
    "                'fix nve all nve',\n",
    "                'fix langevin all langevin %f %f %f %i' % (T, T, Tdamp,\n",
    "                                                           randomseed),\n",
    "                ])\n",
    "\n",
    "    elif integrator == 'nph+l':\n",
    "        start_temp = T*2.+1\n",
    "        Tdamp = 100 * lmp.style.timestep(units)\n",
    "        Pdamp = 1000 * lmp.style.timestep(units)\n",
    "        int_info = '\\n'.join([\n",
    "                'fix nph all nph x %f %f %f &' % (Px, Px, Pdamp),\n",
    "                '                y %f %f %f &' % (Py, Py, Pdamp),\n",
    "                '                z %f %f %f &' % (Pz, Pz, Pdamp),\n",
    "                '                xy %f %f %f &' % (Pxy, Pxy, Pdamp),\n",
    "                '                xz %f %f %f &' % (Pxz, Pxz, Pdamp),\n",
    "                '                yz %f %f %f' % (Pyz, Pyz, Pdamp),\n",
    "                'fix langevin all langevin %f %f %f %i' % (T, T, Tdamp,\n",
    "                                                           randomseed),\n",
    "                ])\n",
    "    else:\n",
    "        raise ValueError('Invalid integrator style')\n",
    "    \n",
    "    return int_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.3. full_relax()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def relax_dynamic(lammps_command, system, potential, mpi_command=None,\n",
    "                  p_xx=0.0, p_yy=0.0, p_zz=0.0, p_xy=0.0, p_xz=0.0, p_yz=0.0,\n",
    "                  temperature=0.0, integrator=None, runsteps=220000,\n",
    "                  thermosteps=100, dumpsteps=None, equilsteps=20000,\n",
    "                  randomseed=None):\n",
    "    \"\"\"\n",
    "    Performs a full dynamic relax on a given system at the given temperature\n",
    "    to the specified pressure state.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    system : atomman.System\n",
    "        The system to perform the calculation on.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    symbols : list of str\n",
    "        The list of element-model symbols for the Potential that correspond to\n",
    "        system's atypes.\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    p_xx : float, optional\n",
    "        The value to relax the x tensile pressure component to (default is\n",
    "        0.0).\n",
    "    p_yy : float, optional\n",
    "        The value to relax the y tensile pressure component to (default is\n",
    "        0.0).\n",
    "    p_zz : float, optional\n",
    "        The value to relax the z tensile pressure component to (default is\n",
    "        0.0).\n",
    "    temperature : float, optional\n",
    "        The temperature to relax at (default is 0.0).\n",
    "    runsteps : int, optional\n",
    "        The number of integration steps to perform (default is 220000).\n",
    "    integrator : str or None, optional\n",
    "        The integration method to use. Options are 'npt', 'nvt', 'nph',\n",
    "        'nve', 'nve+l', 'nph+l'. The +l options use Langevin thermostat.\n",
    "        (Default is None, which will use 'nph+l' for temperature == 0, and\n",
    "        'npt' otherwise.)\n",
    "    thermosteps : int, optional\n",
    "        Thermo values will be reported every this many steps (default is\n",
    "        100).\n",
    "    dumpsteps : int or None, optional\n",
    "        Dump files will be saved every this many steps (default is None,\n",
    "        which sets dumpsteps equal to runsteps).\n",
    "    equilsteps : int, optional\n",
    "        The number of timesteps at the beginning of the simulation to\n",
    "        exclude when computing average values (default is 20000).\n",
    "    randomseed : int or None, optional\n",
    "        Random number seed used by LAMMPS in creating velocities and with\n",
    "        the Langevin thermostat.  (Default is None which will select a\n",
    "        random int between 1 and 900000000.)\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'relaxed_system'** (*float*) - The relaxed system.\n",
    "        - **'E_coh'** (*float*) - The mean measured cohesive energy.\n",
    "        - **'measured_pxx'** (*float*) - The measured x tensile pressure of the\n",
    "          relaxed system.\n",
    "        - **'measured_pyy'** (*float*) - The measured y tensile pressure of the\n",
    "          relaxed system.\n",
    "        - **'measured_pzz'** (*float*) - The measured z tensile pressure of the\n",
    "          relaxed system.\n",
    "        - **'measured_pxy'** (*float*) - The measured xy shear pressure of the\n",
    "          relaxed system.\n",
    "        - **'measured_pxz'** (*float*) - The measured xz shear pressure of the\n",
    "          relaxed system.\n",
    "        - **'measured_pyz'** (*float*) - The measured yz shear pressure of the\n",
    "          relaxed system.\n",
    "        - **'temp'** (*float*) - The mean measured temperature.\n",
    "        - **'E_coh_std'** (*float*) - The standard deviation in the measured\n",
    "          cohesive energy values.\n",
    "        - **'measured_pxx_std'** (*float*) - The standard deviation in the\n",
    "          measured x tensile pressure of the relaxed system.\n",
    "        - **'measured_pyy_std'** (*float*) - The standard deviation in the\n",
    "          measured y tensile pressure of the relaxed system.\n",
    "        - **'measured_pzz_std'** (*float*) - The standard deviation in the\n",
    "          measured z tensile pressure of the relaxed system.\n",
    "        - **'measured_pxy_std'** (*float*) - The standard deviation in the\n",
    "          measured xy shear pressure of the relaxed system.\n",
    "        - **'measured_pxz_std'** (*float*) - The standard deviation in the\n",
    "          measured xz shear pressure of the relaxed system.\n",
    "        - **'measured_pyz_std'** (*float*) - The standard deviation in the\n",
    "          measured yz shear pressure of the relaxed system.\n",
    "        - **'temp_std'** (*float*) - The standard deviation in the measured\n",
    "          temperature values.\n",
    "    \"\"\"\n",
    "    # Build filedict if function was called from iprPy\n",
    "    try:\n",
    "        assert __name__ == pkg_name\n",
    "        calc = iprPy.load_calculation(calculation_style)\n",
    "        filedict = calc.filedict\n",
    "    except:\n",
    "        filedict = {}\n",
    "    \n",
    "    # Get lammps units\n",
    "    lammps_units = lmp.style.unit(potential.units)\n",
    "    \n",
    "    #Get lammps version date\n",
    "    lammps_date = lmp.checkversion(lammps_command)['date']\n",
    "    \n",
    "    # Handle default values\n",
    "    if dumpsteps is None:\n",
    "        dumpsteps = runsteps\n",
    "    \n",
    "    # Define lammps variables\n",
    "    lammps_variables = {}\n",
    "    system_info = system.dump('atom_data', f='init.dat',\n",
    "                              potential=potential,\n",
    "                              return_pair_info=True)\n",
    "    lammps_variables['atomman_system_pair_info'] = system_info\n",
    "    \n",
    "    integ_info = integrator_info(integrator=integrator,\n",
    "                                 p_xx=p_xx, p_yy=p_yy, p_zz=p_zz,\n",
    "                                 p_xy=p_xy, p_xz=p_xz, p_yz=p_yz,\n",
    "                                 temperature=temperature,\n",
    "                                 randomseed=randomseed,\n",
    "                                 units=potential.units)\n",
    "    lammps_variables['integrator_info'] = integ_info\n",
    "    lammps_variables['thermosteps'] = thermosteps\n",
    "    lammps_variables['runsteps'] = runsteps\n",
    "    lammps_variables['dumpsteps'] = dumpsteps\n",
    "    \n",
    "    # Set compute stress/atom based on LAMMPS version\n",
    "    if lammps_date < datetime.date(2014, 2, 12):\n",
    "        lammps_variables['stressterm'] = ''\n",
    "    else:\n",
    "        lammps_variables['stressterm'] = 'NULL'\n",
    "    \n",
    "    # Set dump_keys based on atom_style\n",
    "    if potential.atom_style in ['charge']:\n",
    "        lammps_variables['dump_keys'] = 'id type q xu yu zu c_pe c_ke &\\n'\n",
    "        lammps_variables['dump_keys'] += 'c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]'\n",
    "    else:\n",
    "        lammps_variables['dump_keys'] = 'id type xu yu zu c_pe c_ke &\\n'\n",
    "        lammps_variables['dump_keys'] += 'c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]'\n",
    "\n",
    "    \n",
    "    # Set dump_modify_format based on lammps_date\n",
    "    if lammps_date < datetime.date(2016, 8, 3):\n",
    "        if potential.atom_style in ['charge']:\n",
    "            lammps_variables['dump_modify_format'] = '\"%d %d %.13e %.13e %.13e %.13e %.13e %.13e %.13e %.13e %.13e %.13e %.13e %.13e\"'\n",
    "        else:\n",
    "            lammps_variables['dump_modify_format'] = '\"%d %d %.13e %.13e %.13e %.13e %.13e %.13e %.13e %.13e %.13e %.13e %.13e\"'\n",
    "    else:\n",
    "        lammps_variables['dump_modify_format'] = 'float %.13e'\n",
    "    \n",
    "    # Write lammps input script\n",
    "    template_file = 'full_relax.template'\n",
    "    lammps_script = 'full_relax.in'\n",
    "    template = iprPy.tools.read_calc_file(template_file, filedict)\n",
    "    with open(lammps_script, 'w') as f:\n",
    "        f.write(iprPy.tools.filltemplate(template, lammps_variables, '<', '>'))\n",
    "    \n",
    "    # Run lammps \n",
    "    output = lmp.run(lammps_command, lammps_script, mpi_command)\n",
    "    \n",
    "    # Extract LAMMPS thermo data. \n",
    "    results = {}\n",
    "    thermo = output.simulations[0]['thermo']\n",
    "    \n",
    "    results['dumpfile_initial'] = '0.dump'\n",
    "    results['symbols_initial'] = system.symbols\n",
    "    \n",
    "    # Load relaxed system from dump file\n",
    "    last_dump_file = str(thermo.Step.values[-1])+'.dump'\n",
    "    results['dumpfile_final'] = last_dump_file\n",
    "    system = am.load('atom_dump', last_dump_file, symbols=system.symbols)\n",
    "    results['symbols_final'] = system.symbols\n",
    "    \n",
    "    # Only consider values where Step >= equilsteps\n",
    "    thermo = thermo[thermo.Step >= equilsteps]\n",
    "    results['nsamples'] = len(thermo)\n",
    "    \n",
    "    # Get cohesive energy estimates\n",
    "    natoms = system.natoms\n",
    "    results['E_coh'] = uc.set_in_units(thermo.PotEng.mean() / natoms, lammps_units['energy'])\n",
    "    results['E_coh_std'] = uc.set_in_units(thermo.PotEng.std() / natoms, lammps_units['energy'])\n",
    "\n",
    "    results['E_total'] = uc.set_in_units(thermo.TotEng.mean() / natoms, lammps_units['energy'])\n",
    "    results['E_total_std'] = uc.set_in_units(thermo.TotEng.std() / natoms, lammps_units['energy'])\n",
    "    \n",
    "    results['lx'] = uc.set_in_units(thermo.Lx.mean(), lammps_units['length'])\n",
    "    results['lx_std'] = uc.set_in_units(thermo.Lx.std(), lammps_units['length'])\n",
    "    results['ly'] = uc.set_in_units(thermo.Ly.mean(), lammps_units['length'])\n",
    "    results['ly_std'] = uc.set_in_units(thermo.Ly.std(), lammps_units['length'])\n",
    "    results['lz'] = uc.set_in_units(thermo.Lz.mean(), lammps_units['length'])\n",
    "    results['lz_std'] = uc.set_in_units(thermo.Lz.std(), lammps_units['length'])\n",
    "    results['xy'] = uc.set_in_units(thermo.Xy.mean(), lammps_units['length'])\n",
    "    results['xy_std'] = uc.set_in_units(thermo.Xy.std(), lammps_units['length'])\n",
    "    results['xz'] = uc.set_in_units(thermo.Xz.mean(), lammps_units['length'])\n",
    "    results['xz_std'] = uc.set_in_units(thermo.Xz.std(), lammps_units['length'])\n",
    "    results['yz'] = uc.set_in_units(thermo.Yz.mean(), lammps_units['length'])\n",
    "    results['yz_std'] = uc.set_in_units(thermo.Yz.std(), lammps_units['length'])\n",
    "    \n",
    "    results['measured_pxx'] = uc.set_in_units(thermo.Pxx.mean(), lammps_units['pressure'])\n",
    "    results['measured_pxx_std'] = uc.set_in_units(thermo.Pxx.std(), lammps_units['pressure'])\n",
    "    results['measured_pyy'] = uc.set_in_units(thermo.Pyy.mean(), lammps_units['pressure'])\n",
    "    results['measured_pyy_std'] = uc.set_in_units(thermo.Pyy.std(), lammps_units['pressure'])\n",
    "    results['measured_pzz'] = uc.set_in_units(thermo.Pzz.mean(), lammps_units['pressure'])\n",
    "    results['measured_pzz_std'] = uc.set_in_units(thermo.Pzz.std(), lammps_units['pressure'])\n",
    "    results['measured_pxy'] = uc.set_in_units(thermo.Pxy.mean(), lammps_units['pressure'])\n",
    "    results['measured_pxy_std'] = uc.set_in_units(thermo.Pxy.std(), lammps_units['pressure'])\n",
    "    results['measured_pxz'] = uc.set_in_units(thermo.Pxz.mean(), lammps_units['pressure'])\n",
    "    results['measured_pxz_std'] = uc.set_in_units(thermo.Pxz.std(), lammps_units['pressure'])\n",
    "    results['measured_pyz'] = uc.set_in_units(thermo.Pyz.mean(), lammps_units['pressure'])\n",
    "    results['measured_pyz_std'] = uc.set_in_units(thermo.Pyz.std(), lammps_units['pressure'])\n",
    "    results['temp'] = thermo.Temp.mean()\n",
    "    results['temp_std'] = thermo.Temp.std()\n",
    "    \n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Run calculation function(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict = relax_dynamic(lammps_command, system, potential,\n",
    "                             mpi_command = mpi_command,\n",
    "                             p_xx = pressure_xx,\n",
    "                             p_yy = pressure_yy,\n",
    "                             p_zz = pressure_zz,\n",
    "                             p_xy = pressure_xy,\n",
    "                             p_xz = pressure_xz,\n",
    "                             p_yz = pressure_yz,\n",
    "                             temperature = temperature,\n",
    "                             runsteps = runsteps,\n",
    "                             integrator = integrator,\n",
    "                             thermosteps = thermosteps,\n",
    "                             dumpsteps = dumpsteps,\n",
    "                             equilsteps = equilsteps,\n",
    "                             randomseed = randomseed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['dumpfile_initial', 'symbols_initial', 'dumpfile_final', 'symbols_final', 'nsamples', 'E_coh', 'E_coh_std', 'E_total', 'E_total_std', 'lx', 'lx_std', 'ly', 'ly_std', 'lz', 'lz_std', 'xy', 'xy_std', 'xz', 'xz_std', 'yz', 'yz_std', 'measured_pxx', 'measured_pxx_std', 'measured_pyy', 'measured_pyy_std', 'measured_pzz', 'measured_pzz_std', 'measured_pxy', 'measured_pxy_std', 'measured_pxz', 'measured_pxz_std', 'measured_pyz', 'measured_pyz_std', 'temp', 'temp_std'])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Report results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1. Define units for outputting values\n",
    "\n",
    "- __length_unit__ is the unit of length to display values in.\n",
    "- __energy_unit__ is the unit of energy to display values in.\n",
    "- __pressure_unit__ is the unit of pressure to display values in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "length_unit = 'angstrom'\n",
    "energy_unit = 'eV'\n",
    "pressure_unit = 'GPa'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2. Display number of samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# of samples =  2001\n"
     ]
    }
   ],
   "source": [
    "print('# of samples = ', results_dict['nsamples'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3. Print mean and standard deviations of measured values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Measured thermo data with standard deviation:\n",
      "Ecoh = -4.413022 +-  0.000507 eV\n",
      "lx   = 35.330025 +-  0.025861 angstrom\n",
      "ly   = 35.329559 +-  0.026324 angstrom\n",
      "lz   = 35.330773 +-  0.026982 angstrom\n",
      "xy   = -0.000836 +-  0.029263 angstrom\n",
      "xz   = -0.000028 +-  0.027548 angstrom\n",
      "yz   =  0.000586 +-  0.024387 angstrom\n",
      "T    = 299.963948 +-  3.788839 K\n",
      "Pxx  =  0.000186 +-  0.190085 GPa\n",
      "Pyy  =  0.000212 +-  0.188994 GPa\n",
      "Pzz  =  0.000974 +-  0.190964 GPa\n",
      "Pxy  =  0.000299 +-  0.125513 GPa\n",
      "Pxz  =  0.000353 +-  0.118209 GPa\n",
      "Pyz  =  0.000150 +-  0.103573 GPa\n"
     ]
    }
   ],
   "source": [
    "print('Measured thermo data with standard deviation:')\n",
    "print('Ecoh = %9f +- %9f %s' % (uc.get_in_units(results_dict['E_coh'], energy_unit),\n",
    "                                    uc.get_in_units(results_dict['E_coh_std'], energy_unit),\n",
    "                                    energy_unit))\n",
    "print('lx   = %9f +- %9f %s' % (uc.get_in_units(results_dict['lx'], length_unit),\n",
    "                                    uc.get_in_units(results_dict['lx_std'], length_unit),\n",
    "                                    length_unit))\n",
    "print('ly   = %9f +- %9f %s' % (uc.get_in_units(results_dict['ly'], length_unit),\n",
    "                                    uc.get_in_units(results_dict['ly_std'], length_unit),\n",
    "                                    length_unit))     \n",
    "print('lz   = %9f +- %9f %s' % (uc.get_in_units(results_dict['lz'], length_unit),\n",
    "                                    uc.get_in_units(results_dict['lz_std'], length_unit),\n",
    "                                    length_unit)) \n",
    "print('xy   = %9f +- %9f %s' % (uc.get_in_units(results_dict['xy'], length_unit),\n",
    "                                    uc.get_in_units(results_dict['xy_std'], length_unit),\n",
    "                                    length_unit))\n",
    "print('xz   = %9f +- %9f %s' % (uc.get_in_units(results_dict['xz'], length_unit),\n",
    "                                    uc.get_in_units(results_dict['xz_std'], length_unit),\n",
    "                                    length_unit))     \n",
    "print('yz   = %9f +- %9f %s' % (uc.get_in_units(results_dict['yz'], length_unit),\n",
    "                                    uc.get_in_units(results_dict['yz_std'], length_unit),\n",
    "                                    length_unit)) \n",
    "print('T    = %9f +- %9f %s' % (results_dict['temp'],results_dict['temp_std'], 'K')) \n",
    "print('Pxx  = %9f +- %9f %s' % (uc.get_in_units(results_dict['measured_pxx'], pressure_unit),\n",
    "                                    uc.get_in_units(results_dict['measured_pxx_std'], pressure_unit),\n",
    "                                    pressure_unit)) \n",
    "print('Pyy  = %9f +- %9f %s' % (uc.get_in_units(results_dict['measured_pyy'], pressure_unit),\n",
    "                                    uc.get_in_units(results_dict['measured_pyy_std'], pressure_unit),\n",
    "                                    pressure_unit)) \n",
    "print('Pzz  = %9f +- %9f %s' % (uc.get_in_units(results_dict['measured_pzz'], pressure_unit),\n",
    "                                    uc.get_in_units(results_dict['measured_pzz_std'], pressure_unit),\n",
    "                                    pressure_unit)) \n",
    "print('Pxy  = %9f +- %9f %s' % (uc.get_in_units(results_dict['measured_pxy'], pressure_unit),\n",
    "                                    uc.get_in_units(results_dict['measured_pxy_std'], pressure_unit),\n",
    "                                    pressure_unit)) \n",
    "print('Pxz  = %9f +- %9f %s' % (uc.get_in_units(results_dict['measured_pxz'], pressure_unit),\n",
    "                                    uc.get_in_units(results_dict['measured_pxz_std'], pressure_unit),\n",
    "                                    pressure_unit)) \n",
    "print('Pyz  = %9f +- %9f %s' % (uc.get_in_units(results_dict['measured_pyz'], pressure_unit),\n",
    "                                    uc.get_in_units(results_dict['measured_pyz_std'], pressure_unit),\n",
    "                                    pressure_unit)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.4. Compute lattice constant with standard error\n",
    "\n",
    "**NOTE**: This step makes two assumptions\n",
    "\n",
    "1. The crystal structure is cubic and remains cubic after relaxation.  Check values above to verify this.\n",
    "\n",
    "2. The thermosteps parameter is large enough that the measurements are not correlated. If thermosteps &ge; 100 this is likely a sound assumption."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cubic lattice constant with standard error:\n",
      "a =  3.533012 +-  0.000034 angstrom\n"
     ]
    }
   ],
   "source": [
    "a = results_dict['lx'] / sizemults[0]\n",
    "b = results_dict['ly'] / sizemults[1]\n",
    "c = results_dict['lz'] / sizemults[2]\n",
    "\n",
    "a_std = results_dict['lx_std'] / sizemults[1]\n",
    "b_std = results_dict['ly_std'] / sizemults[2]\n",
    "c_std = results_dict['lz_std'] / sizemults[2]\n",
    "\n",
    "a_mean = (a + b + c) / 3\n",
    "a_combined_std = ((a_std**2 + b_std**2 + c_std**2 \n",
    "                   + (a - a_mean)**2 + (b - a_mean)**2 + (c - a_mean)**2) / 3)**0.5\n",
    "a_standard_error = a_combined_std * (3 * results_dict['nsamples'])**-0.5\n",
    "\n",
    "print('Cubic lattice constant with standard error:')\n",
    "print('a = %9f +- %9f %s' % (uc.get_in_units(a_mean, length_unit),\n",
    "                             uc.get_in_units(a_standard_error, length_unit),\n",
    "                             length_unit))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.5. Load final configuration and show box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [35.319,  0.000,  0.000]\n",
      "bvect =  [-0.027, 35.360,  0.000]\n",
      "cvect =  [ 0.026, -0.004, 35.285]\n",
      "origin = [-0.159, -0.180, -0.143]\n"
     ]
    }
   ],
   "source": [
    "finalsystem = am.load('atom_dump', results_dict['dumpfile_final'],\n",
    "                      symbols=results_dict['symbols_final'])\n",
    "print(finalsystem.box)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
